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NONLINEAR SCHRODINGER EQUATIONS WITH A 
MULTIPLE-WELL POTENTIAL AND A STARK-TYPE 
PERTURBATION 

ANDREA SACCHETTI 


Abstract. A Bose-Einstein condensate (BEC) confined in a one-dimensional 
lattice under the effect of an external homogeneous field is described by the 
Gross-Pitaevskii equation. Here we prove that such an equation can be re¬ 
duced, in the semiclassical limit and in the case of a lattice with a finite num¬ 
ber of wells, to a finite-dimensional discrete nonlinear Schrodinger equation. 
Then, by means of numerical experiments we show that the BEC’s center of 
mass exhibits an oscillating behavior with modulated amplitude; in particu¬ 
lar, we show that the oscillating period actually depends on the shape of the 
initial wavefunction of the condensate as well as on the strength of the non¬ 
linear term. This fact opens a question concerning the validity of a method 
proposed for the determination of the gravitational constant by means of the 
measurement of the oscillating period. 


1. Introduction 


Laser-cooled atoms have drawn a lot of attention as for potential applications to 
interferometry and high-precision measurements, from the determination of gravi¬ 
tational constants to geophysical applications mrnmm, see also PH E2] for a 
recent review. The idea of using ultracold atoms moving in an accelerated optical 
lattice [2 El EH EH [27] has opened the field to multiple applications. In particular, 
by means of the method proposed by Clade et al [|9j, a value for the constant g 
has been measured using ultracold strontium atoms confined in a vertical optical 
lattice mi; such a result has been improved by using a larger number of atoms and 
reducing the initial temperature of the sample m ■ Determination of g has been 
obtained by measuring the period T of the Bloch oscillations of the atoms in the 
vertical optical lattice; recalling that 


2ttIi 
mgb ’ 


( 1 ) 


where m is the mass of the Strontium atom, H is the Planck constant and b is the 
lattice period, then a precise value of the constant g has been obtained by means of 
the experimental measurements of the oscillating period. Since Bloch oscillations 
with period 0 have been predicted by the Bloch Theorem [8] only for a one-body 
particle in a periodic field and under the effect of a Stark potential then it has been 
chosen, in the experiments above, a particular Strontium’s isotope 88 Sr; in fact, 
the scattering length a s of atoms 88 Sr is very small and thus it has been assumed 
by [12 E0] that the effects of the atomic binary interactions are negligible. The 
obtained value for the constant g was consistent with the one obtained by classical 
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gravimeters; but it was affected by a relative uncertainty of order 6 x 10 -6 because 
of a larger scattering in repeated measurements, mainly due to the initial position 
instability of the trap. Such a technique is also proposed to measure surface forces 
[25] , too. 

The critical point of this experimental procedure concerns the validity of the 
Bloch Theorem and the estimate of the effect of the atomic binary interactions on 
the oscillating period of the BEC. In order to discuss this point here we are inspired 
by a realistic model of a one-dimensional cloud of cold atoms in a periodical optical 
lattice under the effect of the gravitational force. The periodic potential has the 
shape 

V per {x) = V 0 sin 2 (k L x) (2) 

where b = l is the period, and Xl = §p The one-dimensional BEC is governed 
by the one-dimensional time-dependent Gross-Pitaevskii equation with a periodic 
potential and a Stark potential 

ihd t il) = H B ip + fxip + 7 | ip\ 2 ip , / = mg, (3) 

where the wavefunction ip(-,t) £ L 2 (R,dx) is normalized to one: 

II^M)I U 2 = ll^o(-)IU 2 = G 

and where 

Hb = + Vper{x) 

is the Bloch operator with periodic potential V per (x). By 7 we denote the effective 
one-dimensional nonlinearity strength. 

It is a well known fact (see §6.1 by [5j) that when the wavefunction ip is prepared 
on the first band of the Bloch operator and if the nonlinear term is absent, i.e. 7 = 0 , 
then the dominant term of the wavefunction ip exhibits a periodic behavior with 
Bloch period T within an interval with amplitude y|j, where Bi is the width of the 
first band and where f £ ffi. is the strength of the external homogeneous field (in 
the case of / = mg then / takes only positive values, obviously). Therefore, for 
times of the order of the Bloch period T we may assume that the motion of the 
BEC occurs in a finite interval. Hence, we can restrict ourselves to the analysis 
of equation Q in a suitable finite interval and then we may assume to consider a 
multiple-well potential Vn(x) (with a fixed number N of wells) and that the Stark 
potential x is replaced by a Stark-type potential Wn(x) due to an homogeneous 
external field which acts only in a bounded region containing the N wells (see Fig. 
[I]). That is, instead of (|3j) we consider, as a model for a BEC in an optical lattice 
under an external homogeneous field, the time-dependent non-linear Schrodinger 
equation (NLS) 

f ied t ip = H N ip + fW N [x)ip + ~/\ip\ 2 ip, H n = -e 2 d 2 x + V N , , 

\ ip{x, 0) = ipo(x) 1 

where e > 0 plays the role of the semiclassical parameter (we prefer to denote here 
the small semiclassical parameter by e instead of the usual notation h because in a 
subsequent section we'll discuss a real physical model where h will assume its fixed, 
physical value; with such a notation it turns out that the Bloch period is given by 
T = jyp;)- We assume that the N wells have all the same shape and we denote by 
b > 0 the distance between the adjacent absolute minima points. 
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Figure 1 . Plot of the multiple-wells potential Vn (full line) and 
of the Stark-type potential tUy (broken line), where N = 11. By 
b > 0 we denote the distance between the adjacent absolute minima 
points. 


The study of the dynamics of the wavefunction ip, solution of 0 , is then achieved 
by means of a discrete nonlinear Schrodinger equation (DNLS). The idea is basically 
simple and it consists in assuming that the wavefunction ip may be written as a 
superposition of vectors ut{x) localized on the i— th cell of the lattice; that is 

N 

^(x,t) ~ y ^cj(t)ui(x) . 

i=i 

Such an approach has been successfully used in the cases of semiclassical NLS with 
multiple-well potentials [Mj or with periodic potentials (see PUCEElcni]), without 
the external field with potential W jv- Eventually, ue(x) may coincide with the 
Wannier function uY (x) associated to the first band of the Bloch operator Hb or 
with the semiclassical single well ground state eigenfunction it| c (x). By means of 
such an approach the unknown functions c^(t) turn out to be the solutions of a 
system of time-dependent equations which dominant terms are given by (here we 
denote ' = jjf) 

iec £ = -X D c e -/3(q+i +q_i) +'y\\u 0 \\'i*\c i \ 2 ct + fMc e , £=1,...,N (5) 

where Ad is the ground state of a single cell potential and where /3 is the hopping 
matrix element between neighboring sites. In fact, the parameter /3 is expected to 
be such that 4/3 is equal to the amplitude B\ of the first band [25]. In ([5]) we’ll fix 
Co = Cn +i = 0. Equation ([5]) represents a discrete nonlinear Schrddinger equation 
(DNLS). 

Our approach is both semiclassical and perturbative. It is semiclassical in 
the sense that it holds true in the semiclassical regime of e small enough; and it 
is perturbative in the sense that the external field / and the nonlinearity power 
strength 7 must be small when e goes to zero (see Hyp. [3] for details). Under these 
conditions we prove the validity of the N- mode approximation 0 with a rigorous 
estimate of the remainder term for times of the order of the Bloch period. Then, we 




























4 


ANDREA SACCHETTI 


numerically solve the ./V-mode approximation and we compute the oscillating 
period taking into account the nonlinear interaction. In fact, the behavior of 
the wavefunction is not simply periodic in time; it turns out that the center of 
mass (xY = shows an oscillating motion with modulated amplitude. The 

oscillating period turns out to be depending on the nonlinearity parameter strength 
7 and we see that it also depends on the distribution of the initial wavefunction if>o- 
In particular, when i/j 0 is a symmetric wavefunction then the oscillating period is 
almost constant for small 7 and it practically coincides with the Bloch period T ; on 
the other hand when ifo is an asymmetrical function the oscillating period actually 
depends on 7. This fact is in contradiction with the Bloch Theorem (which holds 
true when 7 = 0 ), which implies that the Bloch period T does not depend on the 
shape of the initial wavefunction, and it may explain the relatively large uncertainty 
observed by [20] in their experiments, as discussed in the Conclusions. 

The paper is organized as follows. In Section 2 we derive the DNLS ([5]) from 
the NLS Q in the semiclassical limit e —> 0 for times of the order of the Bloch 
period T with a rigorous estimate of the remainder term. In particular: in §2.1 we 
introduce the assumptions and we recall some preparatory results; in § 2.2 we derive 
the DNSL by making use of some ideas previously given by [2J and adapted to the 
case of multiple-well potential with an external Stark-type perturbation. In Section 
3 we consider a realistic experiment and we compute the wavefunction dynamics by 
making use of the DNLS. In particular: in §3.1 we discuss the validity of the N- 
mode approximation for different values of the parameters; in §3.2 we numerically 
compute the wavefunction for times of the order of the Bloch period. In Appendix 
we write the Wannier functions in terms of the Mathieu functions. 

Notation. Let g be a quantity depending on the semiclassical parameter e. In 
the following 

9 = O (e- s °/ £ ) 

means that for any e* > 0 and any p £ (0, S 0 ) there exists C := C p f * such that 
\g\ < Ce- {So ~ p) / e , Ve<E (0,e*). 

Hereafter, by C we denote a generic positive constant independent of e. 

Let N £ N, then by Njv := {1,2,..., } we denote the set of first N positive 

integer numbers. 

By || • \\lp we denote the norm of the Banach space L P (R), by (•,•) we denote 
the scalar product of the Hilbert space L 2 (M). 

2. Derivation of the DNLS © 

2.1. Assumptions and preliminary results. We consider the time-dependent 
non-linear Schrodinger equation ([4]) where Vn is a multiple-well potential and 
W N (x) is a bounded Stark-type potential. In particular we assume that 

Hypothesis 1 . Letv{x) £ Cq°(]R.) be an even (i.e. v(—x) = v(x)) smooth function 
with compact support with a non degenerate minimum value at x = 0: 


v(x) > v min = n( 0 ), Vz £ R, x ± 0 . 
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The multiple-well potential is defined as 

N 

Vn{x) = ^2,v(x- Xi) 
e=i 

for some fixed N > 1, where xt = (£ — b and where b > 0 is such that 

supp v C (-§,+§). 

Hence, by construction the potential Vat ( x ) has exactly N wells with not degen¬ 
erate minima at x = xn, £ £ N^v- 

Remark 1 . We assume thatv{x) is an even function just for argument’s sake. As 
discussed in Remark^this assumption may be removed. Furthermore, we assume 
that v is a smooth function as usual; in fact, a lessere regularity (e.g. C 2 ) would 
be enough. 

Hypothesis 2 . LetWi\r(x) £ C(R) be the monotone not decreasing function defined 
as 

( —L if x < —L 
Wn(x) = < x if x £ [— L, L] 

L if x > L 

for some L > b . 

That is the Stark-type potential Wn is linear in the region containing the wells 
and it is a constant function outside this region (see Fig. [l]). In the “limit” where 
N goes to infinity the potential Vn becomes a periodic potential with period b and 
the external potential Wn becomes the Stark potential x. 


Remark 2 . We restrict ourselves to a multiple-well potential Vn with a finite 
number of wells only for sake of simplicity; one could consider the case of a periodic 
potential by making use of the tools developed by 0 - On the other side, the 
assumption on Wn is not merely for the sake of simplicity; actually, the Stark-type 
potential Wn is a bounded operator while the Stark potential x is not a bounded 
operator and this fact is a source of several technical problems. In fact, in real 
experiments the BEC are trapped in a finite spatial region. 


Hypothesis 3. We assume to be in the semiclassical limit, that is we look for 
the solution of in the limit of e that goes to zero. We assume also that the 
other two parameters 7 and f are small for e small. That is we assume that there 
exists e* > 0 such that 


C e -( s o-p)/* < I/) < Ce s , Ve £ (0,e*), 


for some s > 2, C > 0 and p £ (0, So) independent of e; furthermore, we assume 
also that 



I/I 


< C 


( 6 ) 


for some positive constant C and for any e £ (0, e*). 


The self-adjoint extension of the linear Schrodinger operator formally defined on 
h 2 (R) as 


Hn — + Vn 
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has an almost degenerate ground state with dimension N. More precisely, let A^, 
£ G Nat, be the lowest eigenvalues of H N with associated normalized eigenvectors 
t>i. In particular we have that (see Lemma 2 [25]) 

\i = \ D - 2/3 cos + 0(e°°)e ~ s °^ e , £ € N n , 

where 



is the Agmon distance between two wells and is the ground state of the single 
well operator —e 2 df. x + v, where the single well potential v has been introduced 
by Hyp. 1. The numerical pre-factor /3 is the hopping matrix element between 
neighboring wells, and it is such that 4/3 is asymptotic to the amplitude of the first 
band of the periodic Bloch operator Hg; he. 4/3 ~ B\ := E[ — E\ where E\ and 
E\ are, respectively, the bottom and the top of the first band. Such a numerical 
pre-factor is going to be exponentially small, i.e. 

/3 = 0{e- So/£ ) as e —)• 0 + . 

Remark 3. Hyp. means that, from a practical point of view, the parameter f 
cannot be arbitrarily small, but it has a lower bound of order /3. On the other hand, 
the parameter 7 may be arbitrarily small. 

The associated normalized eigenvectors are given by [251 

N 

v e = Y J ^ J uf + 0{e 00 )e- s ° /£ 

3=i 

where 

n~ . (■„ * ^ 

“« = “« = VFTI“’b wTTj 

and where u* c ( x) is the semiclassical single well ground state eigenfunction localized 
on the j- th cell; by construction and since v(x) is an even function then 

Uj c (x ) = Uq c (x — Xj) and Uq c (x) = Uq c (—x) . (7) 

Now, let n be the projection operator associated with the N eigenvalues i.e. 

N 

n= ") v t 
1 =1 

and let 


n c = 1 - n. 

Let F = n(L 2 (K)) be the ./V-dimensional space spanned by the N eigenvectors vi, 
i e N n . 

Remark 4. Let a{H^) be the spectrum of H^; then it is a well known semiclassical 
result that 

C-h < dist ({A \ {A/}*,) < Ce 
for some positive constant C > 0. Hence, since Hjy is a self-adjoint operator then 
|| Wn - AD] -1 n c || Ai2 ^ i2) < Ce~ 1 
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for some C > 0. 


Remark 5. By |14] it has been proved that there exists a suitable orthonormal base 
U£, i £ Nn, of the space F. The functions are practically localized on the l-t.h 
well. More precisely, they are such that 

i. || U£ — Mf|| L p = O (e -So//e ) for any p £ [ 2 , +oo] and any l € Njv; 

ii. \\u£Uj\\ L i = O (e- So ^-^/ e ) for any j,£ £ Njv/ 

iii. ||u^||iP < Ce~^v , p £ [2,+oo], and ||3 x ^||l 2 < Ce^ 1 / 2 for any £ £ Njv; 

iv. The matrix with elements (uf,H^Uj) can be written as 

((ue, H N Uj )) = Xd^n - PT + D n 


where T is the tridiagonal Toeplix matrix such that 


T, 




0 if | j - £| 7^ 1 

i */ ii - £| = l 


and where the remainder term D jv is a bounded linear operator from £ p (Njv) 
to £ p (Njv) with bound 

II-C , jv||z:(^p(Njv)->-^(Njv)) = O {e ( s o+“)/<^ , p£ [l,+oo], 


for some a > 0 . 


We finally assume that the initial state is prepared on the first N “ground states”. 
That is 


Hypothesis 4. n^o = 0. 

It is well known that under the assumptions above the NLS Q is locally well 
posed, and the conservation of the norm and of the energy mm 

£(V0 = (V> 5 H N if) + ^HV’llli + f(if, w N %p) 

easily follow: 

IIV’(-,*)I|l 2 = II^MOIIl 2 and £ {${■, t)) = £ (^o(-)) ■ 

Furthermore the following a priori estimate follows, too. 


Lemma 1. There exists a positive constant C > 0 such that 

HV’llff 1 < Ce~ 1/2 and \\if>\\ p LP < C<T ^ , \j v g [ 2 , +oo]. 

Proof. Indeed, from Theorem 2 by [23J and its remarks it follows that 

||V^||i 2 < eVA and II^IIlp < CA s w 
for some C > 0 and e small enough, where 

* _ £(l/>0) Vniin 

~ ? 

and where I X min = min a; [Viv(^) + fWjy(x)]. In particular, since fW^{x) > 
0(e s ) for some s > 2, because L is hxed, and since II c ^o = 0 then A 
therefore 

IIV^IIl 2 < CeT 1 ^ 2 and ||^||lp < Ce~ E ^~ . 


-fL = 


c- 1 . 


□ 


Hence, the global well-posedness of the NLS follows mm- 
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2.2. N-mode approximation. Let ip be the normalized solution of the NLS equa¬ 
tion written in the formula 


N 


ip = ipi+ip c , ipi = Ihp = '^2 c e u i an d V’c = n c ip , 


( 8 ) 


for some complex-valued functions C({t). By substituting § into the NLS 0 then 
it takes the formula 


ied t ct = {ut,H N ipx) +7 {u e , \ip(-, r)\ 2 ip{-, r)) + f(u e ,W N ip(-,T )) 
ied t ip c = H N ip c + "/U c \ip(-,T)\' 2 ip(-,T) + fU c W N ip(-,T) 


(9) 


We are going now to get an a priori estimate of the remainder term ip c . First of 
all we rescale the time t —> r = jt and we redefine the wavefunction up to a gauge 
factor ip(x,t ) —► ip(x,r) := e~ lXl>t I e ip{x,f). The Bloch period becomes 

_ Prr_ W 

TB e I/I 6 

Hence, 0 becomes (where ' denotes the derivative with respect to r) 

i/3^ = ( h n - A D )ipi) +7 {ue, t)\ 2 ip{-, t)) + }{u t ,W N ip(-,t)) . . 

iplp'c = (H N - \ D )lP c + 7 U c \lP(;t)\ 2 lP(;t) + fn c W N lP(;t) 1 UJ 

Theorem 1. Let Hyp.l-f be satisfied; then it follows that the remainder ip c can be 
estimated for times of order of the Bloch period. That is for any fixed M £ N it 
follows that 

max \\iP c (-,t)\\ L 2 < C— 
tg[o,Mt b ] e 

for some positive constant C > 0. 


Proof. From the first equation of (10) and recalling that 
N 


M r )i 2 = n^iiii 2 = 1 - w^wh < i 


then a priori estimate follows 


Wt\ < ; Wl ? + h!wii. 


i/i 

p 


IM; 


( 11 ) 


because ||m^||l 2 = 1 and ||V’||l 2 = 1, and from Remark[5]iv. and Lemmajl] 
Concerning ip c it satisfies to the following integral equation 


where we set 


ip c = I + II 

J e -i(i? N -A D )(T- 5 )//3 H c Ads 
II := -- [ e -i(^iv-A o){r-s)/P n Bds 

P Jo 


1 = “j 
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and where A and B are defined as 
A := 7|^i| 2 i/'i + fW N ipi 

B := 7 [V>iV>c + |V’c| 2 V’c + 21-01 |Vc + 2 |V’ c | 2 V’i + V’lV’c] + fW N ip c 

such that 

A + B = 'y\’ip\ 2 ip + . 

By means of standard arguments P2 and making use of the fact that the operator 
Wn is bounded then it follows that 

Lemma 2. Let 


r = M 


,- 1/2 . 


1 / 1 - 


Then the functions A and B are such that 

\\A \\ L 2 < CT, \\B \\ L 2 < CT||V» c ||i 2 and 


dA 


dr 


< cr 2 /r 1 . 


L 2 


Proof. Indeed, 

\\A\\v < M HV'illiooll^ilU- + 1/111^11^11^11^ < C 


."I / 2 


I/I 


since ||'/’i||l“ < Crnax^ ||m^||l“ < Ce 1 ' 4 ‘. Similarly, the estimate of the function 
B follows recalling that ||V’c||l“ < Ce -1 / 4 from Lemma [l] Finally, the estimate 


concerning JA immediately follows from ( 11 ). 

Hence, the estimates of the integrals I and II follow; in particular, integral II can 
be simply estimated as 


□ 


\\II\\ L 2 < err 1 flKMIMa 
Jo 


since 


,—i(H N -\ D )(T-s)/p 


= 1 . 


C ( L 2 ^ L 2 ) 

On the other hand, before to get the estimate of integral I we perform an integration 
by parts in order to gain a pre-factor /3: 


I = 


—j e —A dHt-sV/3^ _ x d ]-^u c A ° + 

+i f _ XD^Uc^ds 

Jo 9s 


From this fact and recalling that (Remark 4) 

II 1 n c ||£(L 2 — s-L 2 ) < Ce 1 

then 


| [2 < Ce max 

«e[0,r] 


■ 

dA 

~ 

II^IU 2 +T 

ds 

L 2 - 


< Oe- 1 r[l+r/3- 1 r]. 


Therefore, we have that 

||^ciu= < crp- 1 [ T \\M; s )|| L 2 ds + Ce- 1 r( 1 + rrV). 
Jo 

From the Gronwall’s Lemma it follows that 

\M; t)\\ L 2 < C'e- 1 r(l + r p- l T)e CT P- lT 
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In particular we observe that 

max U(-,t)\\ L 2 < Ce~ 1 T(l + /3- 1 TMT B )e cr P- lMTB < CTe" 1 < C ^1 
re[o,Mr B ] e 


proving the Theorem since T < C\f\ from Hyp. 3 and since t b = — — 


b I/I' 


□ 


We are going now to estimate the solutions ce of the first equation of (101 which 
can be written as 


i/3c' e = (ue, (H n - \ D )^i) + (ue, A) + (ue, B) 

where the term (ue, (Hn — Ad)^i) can be represented by property iv. of Remark[5] 
Concerning the term (ue, B) the following estimate uniformly holds with respect to 
the index £ 


\(ue,B)\ < \\B\\ La < CT\W c \\ l3 . 

Furthermore 

N N 

(ue,A) = 7 CjC k Cm(ue,UjU k u m ) + f^2cj(ue,W N Uj) 

j,k,m =1 j =1 

= 7 |q \ 2 ceWueW^i + fc e (ue, W N ue ) + 7 r“ + fr\ 

where 

r a t = 22 cjc k c m (ue,uju k u m ) 

: | j—£\-\-\m—£\ + \k— ^|>0 

and 

r b e= 22 Cj(ue,W N Uj) 
j/eN jv , j^t 

are remainder terms. 

We have that 


Lemma 3. The following estimates uniformly hold with respect to the indexes 
l, j , m and k: 


i. (ue, W N ue ) = tb + 6 (e _s °/ e ); 

ii. (ue, W n u 3 ) = 6 (e-^b’-^IA); 

iii. (ue, UjU m Uk) = O (e~ s ° r / e ) where 


r = max [\j - £\, \m - £\, |k - £\, \j - m\, \j - k\, \k - m\]. 


Proof. Indeed, let Ie = [xe — b, Xe + b], then 

(ue,W N ue) = / \ue(x)\ 2 xdx + / \ue(x)\ 2 W N (x)dx 

Jlf JR\If 


where ue(x) = u 0 (x — xf) + O (e s ’°/ e ) from (M) and Remark 4. Therefore 


>h 


\ue(x)\ 2 xdx = £b \u 0 (x)\ 2 dx + j \u 0 (x)\ 2 xdx + O 
Jlo 


= £b 


Ji 0 

l U o|||2 — ll U o|||2(R\/ 0 ) 


( e -5o/,) 


+ 0 


( e -SoA) 


where uq is normalized and fj \uq(x )\ 2 xdx = O (e s °/ e ) because uq(x) = uq(—x) + 
O (e -50 / 6 ). From this fact and since ||||z, 2 (r\z £ ) = 0(e~ s °^ h ) ( see Lemma 4 iii. 
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and Lemma 5 by fl4j) then the asymptotic behavior i. follows. The other two 
asymptotic behaviors ii. and iii. similarly follow from property ii. by Remark [5] 
indeed 

\(ue, W N Uj)\ < HWjvlU-.Hu/UillLi = 6 
and, where we assume that r = \j — t\, 

\(ue,uju m u k )\ < ||u m |U~||u||L°o||u/Ui|Ui = 6 (e~ s ° r/e ^ 
proving so the estimates ii. and iii.. □ 


From this Lemma and from the previous computation it follows that the first 
equation of (10) becomes a DNLS of the form 


N 

i/3c' e = ^ Ti, m c m + 7 IMII 4 \ce\ 2 cz + fibci + 6 (e~ So/e ) (12) 

m =1 


where 



= \\uo\\ 4 Li +d(e- So/t ) 


and where the remainder terms O (e ' s °/ £ ) are uniform with respect to the index 

L 


Remark 6. In fact, if v(x) is not an even function then by means of standard 
semiclassical arguments it follows that property Lemma 3 i. becomes 


(ue, Wnw) = £b + c+ O (e s °/ e ^ 


for some constant c independent of the index I. In such a case we must add the 
term fccg to the right hand side of the DNLS above and, by means of a gauge choice 

■ f c f 

ci —> cee~ l p T , we can remove this term obtaining again equation (12). 

Now, we are able to prove that 
Theorem 2. Let di{r) be the solutions of the DNLS 

N 


i(3d' e = -pYl Te ’ mdm + \di\ 2 di + flbdi 


(13) 


m =1 


satisfying to the initial conditions di(0 ) = q(0 ), where q(t) and if c are the solutions 
of \10 1). Then, for any fixed M eN it follows that 


max \ci(t) - dt{r)\ = O 

t£[0,Mt b ],1=1,2,...,N 


(e- s °/ £ ) 


0 . 


Proof. The proof is a simply consequence of equation (12) and from the fact that 

tb = xf and H yp- 3 - 


□ 
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3. Numerical analysis of a real model 


We consider the experiment where a cloud of ultracold Strontium atoms 88 SY 
are trapped in a one-dimensional optical lattice with potential ([2]). Realistic data 
for the experiment are [ 20 ] : 

- Lattice period: b = Al/2 = 266 nm, Xl = 532 nm; 

- Lattice potential depth: Vo = A 0 • Er where Er is the photon recoil energy 
Er = 2 ^ A 2 = 50.38 kHz ■ h and where A 0 is between 3 and 10; 

- Mass of the strontium 88 isotope: m = 87.91 au = 1.46 ■ 10 -22 gr ; 

- Effective one-dinrensional nonlinearity strength: let 73 d = 4A/ be the 
effective nonlinearity strength for the three-dimensional Gross-Pitaevskii 
equation, then it is expected that the effective one-dimensional nonlinearity 
strength 7 is of the order :26j 


73 D 
7 ~ 2ird 2 ± 

where d± is the oscillator length of the transverse confinement; here a s 
denotes the scattering length of the Strontium 88 isotope: a s = —ao7-13oo, 
where ao is the Bohr radius; J\f is the number of atoms of the condensate; 
in typical experiments d± ss 180 • 10 -6 m and A f = 10 5 -7 10 6 ; 

- Acceleration constant g = 9.807to/s 2 . 


The confined BEC is governed by Eq. Q and here we make use of the IV-mode 
approximation (13), that is the wavefunction ij) has the form if} ~ Ye where ce 
are the solutions of (13) and where Uf are functions localized on the £-th lattice site. 
In order to justify the validity of such an approximation we’ll check if the model 
is in the semiclassical regime , that is if the first band is almost flat and if semi- 
classical approximation u| c agrees or not with the Wannier function ■ Such a 
qualitative criterion has been also adopted by other authors mmm and we’ll see 
that our results agree with the results contained in these papers. In particular, in 
mi has been computed the hopping matrix elements (iq, H^Uj) too, where it has 
been numerically verified that these coefficients are negligible when | j — £| > 1 for 
A 0 > 10; thus, for such values of A 0 it is admitted that the IV-mode approximation, 
consisting to describe Q in terms of a nearest-neighbor model (13), works. 


3.1. Validity of the semiclassical approximation. The semiclassical approxi¬ 
mation Uq c (x) of the wavefunction has dominant behavior 


““ (I > = 


(14) 


in the semiclassical limit, where g = d = 2Vo, Vo = A 0 Er; it is normal¬ 

ized |K c || L 2 = 1 . We may remark that the effective semiclassical parameter in 
adimensional units is given by 

1 2ir 2 h 

7T 0 = b 2 ^Tfi' 


and then the semiclassical approximation may be written as 


c (z) = 


'27t v /A 0 ' |1/4 
b 2 


-x 2 K 2 VXo/b 2 





















13 


Hence 


4 

L 4 


m\x 

1/4 

[^ tiAo 74 

(nh) 2 

\ 

/ 2 b 


We’ll see that for Ao “large enough” (i.e. Ao > 10) then the first band is almost 
flat and the semiclassical function u(f well approximates the Wannier function , 
as we expect to observe in the semiclassical limit A 0 —>• oo (see, e.g., [30]). 


Remark 7. By the scaling 

x—> 2 Izlx , t —> Eftt/h , if(x) 

and setting 


F = 


mg 

2 E R kL 


c = 


7 


2 ErIcl 


y/2k£ 


e = 






then 0 takes the form 


id t if = ~d 2 xx if + sin 2 (a:/2) + Fxif + CIV'lV > IIV’II^ = 1 ■ 


(15) 


Equation (15) is equivalent, up to a change of scale of the time, to the equation 
ied t 4> = -e 2 d 2 x ip + sin 2 if + fxif + j\if\ 2 if 


where we set 


t —> t/\fe , / — Fe 2 , 7 — e 2 ( ■ 


_ 1/2 

and where e = A 0 ' plays the role of the semiclassical parameter. 

We compute now the band functions and the Wannier functions for different 
values of Ao. The semiclassical wavefunction uff is computed by (14), while the 
Wannier function u(x) may be computed by means of the Mathieu functions (see 
Appendix). 


3.1.1. Model A 0 = 3. The first bands of the Bloch operator Fl B = ~^d xx + 
AoE B sm 2 (kLx) have endpoints 
n=l) E\ = 1.43 • E r and E\ = 2.11 • E R ; 

n=2) El = 2.86 • E R and E\ = 5.49 • E R - 

n=3) E\ = 5.56 • E R and E 3 = 10.51 • E R . 

Hence, the values of the width of the first two bands are given by 

Hi := E\-E\ = 0.68 • E R and B 2 := E\-E{ = 2.63 • E R . 

Furthermore it follows that the first gap has amplitude g\ = E% — E\ = 0.77 • E R 
of the order of the first band amplitude, while the width of the other gaps are very 
small (see Fig. [2j left hand side panel). If we compare the first Wannier function 
uft' (x) and the semiclassical approximation u(f(x) it turns out that (see also Fig. 
[ 2 j right hand side panel) 

IIC -<\\h= 0-091 
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Figure 2. Here we plot in the left hand side panel the first three 
band functions E n (k), n = 1, 2,3 and k € ^, +f ], for the Bloch 

operator Hg where Ao = 3. It turns out that the width of the 
first gap is of the same order of the width of the first band. In 
the right hand side panel we plot the graph of the functions u 3 c 
(broken line) and iiff (full line). 


3.1.2. Model Ao = 10. The first bands of the Bloch operator Hg have endpoints 
n=l) E\ = 4.32 • Eg and E{ = 4.58 • Eg; 

n=2) E\ = 7.02 • Eg and E\ = 8.87 • Eg; 

n=3) E\ = 9.54 • Eg and E f 3 = 14.07 ■ Eg. 

Hence, the values of the width of the first two bands are given by 

Bi := E\ — E\ = 0.26 ■ Eg and B 2 := E b 2 - E{ = 1.85 • Eg . 

Furthermore it also follows that the first gap has amplitude g\ = E^ — El = 2.44 -Eg 
is much larger than the amplitude of the first band and that the width of the other 
gaps are very small (see Fig. [3j left hand side panel). If we compare the first 
Wannier function vff (x) and the semiclassical approximation u s 0 c (x) it turns out 
that (see also Fig. [3| right hand side panel) 

114^ -VoWb = 0.055. 

Hence, we may conclude that for Ao = 10 the iV-mode approximation properly 
works. 


3.2. Numerical analysis of the model for A 0 = 10. We have seen that for 
Aq > 10 the IV-mode approximation is justified. For Aq = 10 we have that 


P 


;Bl 


Equation (13) takes the form 


0.065 • Eg. 


N 

id' e = -^2 Te, m d m + v\de\ 2 d( + £Sd {,, l € Nat , 

m—1 
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Figure 3. Here we plot in the left hand side panel the first three 
band functions E n (k), n = 1, 2,3 and k € [—|-, +f ], for the Bloch 
operator Hg where Ao = 10. It turns out that the first band is 
almost flat, in fact its width is 1 / 10 -th of the width of the first gap. 
In the right hand side panel we plot the graph of the functions Uq c 
(broken line) and uff (full line). 


where we set 

7 ||ug c ||^ 4 4J\fTra s H 2 ttAI /4: 1 1 


V = 


P 


m 


b 2 'kcP, ft 


- = -0.151 • 10- 


and 


S=A = ^b =1 . im , 


P P 


The Bloch period is given by 


0.197 


2irh 

T = -- = 1.740 ms. 

mgb 

Hence, the parameters /, 7 and /3 are in a suitable range as discussed in Remark [3] 
Furthermore, the motion of the Bloch oscillator occurs in an interval with width 


Bi 

I/I 


0.26 • Eg m—7 ^ o r u 

-= 9.65 -10 m ~ 3.6 • 0 . 


(16) 


Hence, the ./V-mode approximation with = 40 properly works. 

We consider three different situations. In the first one we assume that the state 
is initially prepared on a single lattice site, that is ipo is a Wannier type function. 
In the other two cases we assume that the initial wavefunction ipo is a symmetric 
or asymmetrical wavefunction initially prepared on different lattice sites. 


3.2.1. ipQ is initially prepared on a single lattice cell. We consider a numerical ex¬ 
periment where tpo{x) = uq(x), that is c^(0) = 0, for l ^ IV/2, and c/v/ 2 (0) = 1 
(where N = 40). In fact, in such a case we observe a breathing motion for the 
wavefunction; that is, the wavefunction, initially prepared in a Wannier state lo¬ 
calized on a single site of the optical lattice, symmetrically spreads in space and it 
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periodically returns to its initial shape (Fig. [4j top panel, obtained for 77 = 0.2). 
Then the expected value of the center of mass 

(x)‘ = 

is practical constant (a;)* « 0 up to small fluctuations. 



t 

Figure 4. In the top panel we plot the absolute value of the 
wavefunction ip(x,t) initially prepared on a single Wannier state 
for 77 = 0 . 2 , it turns out that it symmetrically spreads in space 
and periodically returns to its initial shape without motion of the 
center of mass. In the bottom panel we plot the absolute value 
of the wavefunction initially prepared on several lattice sites for 
77 = 0 . 2 ; it turns out that the center of mass oscillates with no 
marked changes of the shape of the wavefunction. Here T denotes 
the Bloch period and b is the distance between two adjacent wells. 
Dark regions mean that \ip(x,t)\ is practically zero there, white 
regions mean that t)| has its maximum value there. 
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0 

II 

0 

cio = 0.396 ■ 10~ a 

c 2 i = 0.429 

c 3 i = 0.898 • 10 ~ 4 

ci = 0 

c n = 0.151 • 10 ~ 2 

c 22 = 0.347 

c 32 = 0.177 • 10 ~ 4 

S> 

II 

0 

C 12 = 0.502 • 10 ” 2 

C 23 = 0.244 

c 33 = 0.303 • 10 -b 

c 3 = 0 

C 13 = 0.149 • 10 ” 1 

c 24 = 0.149 

C 34 = 0 

O 

II 

cm = 0.363 ■ 10 ” 1 

c 25 = 0.788 • 10 1 

C 35 = 0 

c 5 = 0 

C 15 = 0.788 • 10 ” 1 

c 26 = 0.363 ■ 10 1 

C 36 = 0 

c 6 = 0 

ci 6 = 0.149 

c 27 = 0.149-lO" 1 

C 37 = 0 

c 7 = 0.303 • 10 “ 5 

C 17 = 0.244 

c 28 = 0.502 • 10 -2 

C 38 = 0 

c 8 = 0.177- 10 " 4 

cis = 0.347 

c 29 = 0.151 • 10 -2 

C 39 = 0 

c 9 = 0.898 • 10~ 4 

C 19 = 0.429 

c 30 = 0.396 • 10" 3 

0 

0 

II 

O 


C 20 = 0.460 




TABLE 1. Initial values of the coefficients Cf := q( 0) of the wavefunc- 
tion. The initial wavefunction 0o has a symmetric shape and its width 
is of order of several lattice periods. 


3.2.2. t/jg is a symmetric wavefunction initially prepared on different lattice cells. 
We consider a numerical experiment where N = 40 and tf>o(x) = Y^e= oQ^Wi 
where Cf, have a symmetric Gaussian-type distribution around £ = N/2. That is 
the initial value of the coefficients C((t) are given in Tablejlj the initial wavefunction 
■00 is plotted in Fig. [5j left hand side panel. In such a case the center of mass 
(xf oscillates in space and the wavefunction moves with no marked changes in 
shape (see Fig. [4j bottom panel). In particular, the function (xf exhibits, for 
?7 -=f 0 , an oscillating motion where the wavefunction amplitude is modulated (see 
Fig. [§ and where the oscillating (pseudo-)period (that is the time interval between 
two consecutive minima or maxima points) depends on 77 . In Fig. [7]we plot the 
mean value of the oscillating period of the motion of the center of mass after 14 
oscillations for 77 in the range [— 0 . 1 , + 0 . 2 ]; it turns out that the relative uncertainty 
with respect to the Bloch period is of order 2.4 • 10 -5 . 


3.2.3. 0o Is an asymmetrical wavefunction initially prepared on different lattice 
cells. We consider a numerical experiment where N = 40 and 0o(^) = c z u i( x )i 
where ca have an asymmetrical Gaussian-type distribution. That is the initial value 
of the coefficients c^(t) are given in Table [2j the initial wavefunction is plotted in 
Fig. § right hand side panel. As in the symmetric case the center of mass (xf 
oscillates in space and the wavefunction moves with no marked changes in shape. 
Even in such a case the function (xf exhibits, for 77 0, an oscillating motion where 

the wavefunction amplitude is modulated. In contrast with the symmetric case 
the oscillating (pseudo-)period (that is the time interval between two consecutive 
minima or maxima points) actually depends on 77 ; in Fig. [7] we plot the mean value 
of the oscillating period of the center of mass after 14 oscillations for 77 in the range 
[— 0 . 1 ,+ 0 . 2 ] and it is not almost constant like in the previous case, in particular it 
turns out that the relative uncertainty with respect to the Bloch period is of order 
4.6 • 10~ 4 , which is 20 times the relative uncertainty observed in the symmetrical 


case. 
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Figure 5. Here we plot the absolute value of the initial wave- 
function ijj o prepared on several lattice sites; the left hand side 
panel corresponds to the symmetric initial wavefunction, the right 
hand side panel corresponds to the asymmetrical one. 


Co = 0 

do = 0.180 • 10 ” a 

c 2 i = 0.252 

c 31 = 0.175 • 10 ~ 4 

Ci = 0 

cn = 0.814 -10~ a 

c 22 = 0.170 

c 32 = 0.323 • 10~ b 

ft 

to 

II 

O 

ci 2 = 0.330 • 10 ~ 2 

c 23 = 0.103 

c 33 = 0 

c 3 = 0 

ci 3 = 0.121 • 10" 1 

c 24 = 0.546 • 10 " 1 

C 3 4 = 0 

o 

II 

o 

ci 4 = 0.414 • 10 ” 1 

c 25 = 0.257 • 10 " 1 

c 3 5 = 0 

c 5 = 0 

Ci 5 = 0.133 

c 26 = 0.106 • 10 1 

C 3 6 = 0 

c 6 = 0 

cio = 0.351 

c 27 = 0.386 • 10 -2 

C 3 7 = 0 

C 7 = 0 

C 17 = 0.496 

c 28 = 0.123 • 10 —2 

C 3 8 = 0 

c 8 = 0.614 • 10 ~ 5 

cis = 0.471 

c 29 = 0.340 • 10 a 

C 3 9 = 0 

c 9 = 0.354 • 10 ~ 4 

ci 9 = 0.411 

c 30 = 0.826 • 10 -4 

ft 

O 

II 

o 


C 20 = 0.336 




TABLE 2. Initial values of the coefficients a := Cf(0) of the wavefunc¬ 
tion. The initial wavefunction i/jo has an asymmetrical shape and its 
width is of order of several lattice periods. 


4. Conclusion 


In this paper we have proved that in the semiclassical limit the TV-mode approx¬ 


imation (13), corresponding to a discrete nonlinear Schrodinger equation with a 
finite number of modes, gives the solution of the Gross-Pitaevskii equation Q for a 
BEC in a multiple-well lattice in a Stark-type external field. Furthermore, we have 
numerically solved the TV-mode approximation considering a real model, where for 
some values of the physical parameters the validity of the iV-mode approximation 


(131 seems to be justified. In particular, we have seen that a state initially prepared 


on several wells have an oscillating behavior with modulated amplitude, the oscil¬ 
lating (pseudo-)period is computed for different values of the nonlinear strength 
and it turns out that such a period is practically constant when the initial state 
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Figure 6. Here we plot the motion of the center of mass of the 
wavefunction initially prepared on several lattice sites. The initial 
wavefunction is a symmetric function. Top panel corresponds 
to the case of 77 = 0 . 1 , bottom panel corresponds to the case of 
77 = 0.2. The center of mass rapidly oscillates with modulation of 
the amplitude. The width of the escillations is in a range lesser 
or equal to 36 is agreement with (16). 


is a symmetric one; on the other side, such a period actually depends on the non¬ 
linear strength when the initial state is an asymmetrical one. This observation 
opens a question about the validity of the method proposed by Clade et al jjfj for 
the deterimantion of the gravitational constant g by means of the measurement of 
the oscillating period [121 ED], where it has been assumed that the oscillating pe¬ 
riod coincides with the Bloch period T independently from the shape of the initial 
wavefunction and of the value of the nonlinearity strength parameter. 
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Figure 7. Here we plot the mean value of the pseudo-period of 
the oscillating motion of the center of mass after 14 oscillations, 
as function of the effective nonlinearity parameter 77 . Broken line 
corresponds to the case of a symmetric wavefunction prepared on 
several lattice sites; it turns out that in such a case the oscillating 
period is almost constant. Full line corresponds to the case of 
an asymmetrical wavefunction prepared on several lattice sites; it 
turns out that it actually depends on 77 . Here T denotes the Bloch 
period, while t denotes the oscillating period. 


Appendix A. Band Functions and Wannier functions 


For a generic one-dimensional Bloch operator Hb the spectrum is given by a 
sequence of infinitely many closed intervals named bands. These intervals are the 
image of functions named band functions. The band functions of Hb are denoted 
by E n (k), where the quasimomentum k runs in the Brillouin zone [— ^, + ?]. The 
spectrum of the Bloch operator Hb is given by the bands a{Hs) = U^ =1 [E^, E^] 
where 


Et = 


E n ( 0) if n is even t _ / E n (n/b) if n is even 

E n (ir/b ) if n is odd ’ n ( E n { 0) if n is odd 

In the case of potential ([ 2 ]) the band functions may be explicitly computed. In 
particular let us look for the Bloch functions of the equation 


If we set 


h 2 d 2 

H b 4> = Eip, H b = - — ^2 + V 0 sin 2 {k L x). 


1 , 


P _ 1 - T/ | 2m - V 0 to _ mA 0 E r _ 1 2 

8 [ E 2 Vo ) h 2 ’ q 2klj! V ° h 2 h 2 2 A ° kL 

and recalling that sin 2 (0) = | [1 — cos(20)] then the Mathieu equation takes the 
form 

d 2 


Hn-E 


if> = 0 where H B = — — Vo cos(ga;). 

dx 2 


(17) 
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Figure 8 . Here we plot the motion of the center of mass of the 
wavefunction initially prepared on several lattice sites. The initial 
wavefunction is a symmetric function. Top panel corresponds 
to the case of 77 = 0 . 1 , bottom panel corresponds to the case of 
77 = 0.2. The center of mass rapidly oscillates with modulation 
of the amplitude. The width of oscillation is in a range lesser or 


equal to 36 is agreement with (161. 


It has a fundamental set of solutions [I] 


Mx,£) = c 


AS V 0 1 

-Y, - 2 ’ o qx 

q z q- 2 


and ip 2 (x, S) = -S 

q 


AS V 0 1 


qx 


where S and C denotes the two Mathieu’s functions, satisfying the conditions 


, fn C \ 1 ^1(0, £) n A , (n 3 t/> 2 ( 0 ,£:) 

V’i(0,t) = l,-—-=0 and t/> 2 (0,£)=0,-—-= 1. 
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Hence, the band functions £„ (, k ) associated to the spectral problem © are the 
solutions of the equation /i (£) = cos (kb) where 


V(£) = ipi(b,£) = C 


4£ Vo 
q 2 ’ q 2 ’ n 


Let A = e lkb , then the equation p(£) = cos (kb) can be written as p(£) = 
\(A + A -1 ). We observe that for k £ [0, f] then sin(fcfe) = y/l — p 2 [£). The 
Bloch function is given by m 

X(x,£) 


ip(x,£) = 


VW) 


where 


and 


x(x,£) = il> 2 (b,£)i>i(x,£) + -[A(£) - A (, £)]ip 2 {x,£) 

= i/j 2 (b,£)4>i(x,£) + iy/l - n{£) 2 i> 2 (x,£) 


We recall that the Bloch function ip n (x,k) = ip(x,£ n (k)), where £ n is the band 
function associated to Hg , is normalized to one: 

rb 

1 1pn(x, k)\ 2 d,X = 1 


2n 

T 


and furthermore it is such that 


if) n (x,-k) = ip n (x,k). 

Finally, the Wannier function on the zero-tli cell associated to the n-th band is 
given by 


( b A 1/2 f + ^ /h ( b \ 1/z f 

(x) = (-) J M*,k)dk = 2 (-j l 


b \ 1//2 r+ n / b 


\Ripn ^x , k)dk 


= 2 


b A 1/2 r /b ip 2 (b,£n{k))ipi(x,£ n (k)) 


2t t) 


o y/N(£ n (k)) 

b r /b \/^ 2 {b,£ n (k))'ijj\(x,£ n {k)) 


dk 


\[2'k Jo 
1 


d^i(£ n (k)) 

d£ 


dk 


£^/») ^Mbj£)J>i(x,£)ypW 


d£ 


Js„( 0) y/l - 

since the Mathieu functions are real valued when their arguments are real numbers. 
In particular, 


o 0) : = wi(x) = 


1 


r £ l \JMb,£)ipi(x,£)\f-^jp- 


d£. 


Jei ^l-p 2 (£) 
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